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1.  Introduction 


The  application  of  photographic,  CCD,  and  other  forms  of  imaging  for  the 
purpose  of  estimating  flow  velocities,  has  been  investigated  by  many  researchers  in 
fields  ranging  from  fluid  mechanics  to  vision  research.  In  the  most  common  methods 
for  measiuring  fluid  flow  velocities,  the  flow  is  seeded  with  pai'ticles,  or  markers,  that 
can  be  easily  imaged  and  tracked.  An  extensive  review  of  methods  using  particle 
and  speckle  images  for  fluid  flow  measurement  is  presented  by  Adrian  (1991).  The 
estimation  of  the  motion  and  deformation  of  solids  is  closely  related  to  that  of  fluids. 
A  method  of  determining  displacements  and  stress  intensity  factors  in  solids,  using 
white  light  speckle  images  and  image  correlation  techniques,  is  presented  in  McNeill 
et  al.  1987.  In  the  absence  of  particles,  flows  have  also  been  tagged  with  a  line 
or  grid,  e.g.,  using  laser-induced  photochemical  reactions  (Falco  &  Chu  1987),  or 
laser-induced  fluorescence  (Miles  et  al.  1989).  When  this  is  not  possible,  one  can 
use  markers  that  occur  naturally  in  the  flow,  e.g.,  Bindschadler  and  Scambos  (1991) 
have  correlated  the  translation  of  distinct  surface  features  in  ice  flows  to  determine 
flow  velocities. 

Determining  motion  from  successive  images  is  £dso  of  interest  in  ainimation,  as 
well  as  the  study  of  biologiccd  and  robotic  vision.  Most  investigations  along  these 
lines  have  taken  the  form  of  extracting  the  motion  of  objects  in  an  image  aind,  as  a 
consequence,  they  focus  on  the  motion  of  rigid  objects  and  their  representations.  See 
Hildreth  k.  Koch  (1987)  for  a  review  and  Murray  Sc  Buxton  (1990).  This  approach 
is  somewhat  different  from  the  interests  of  Fluid  Mechanics  where  the  object  of 
interest  is  a  fluid,  highly  deformable  and  often  compressible.  Nevertheless,  mainy 
results  from  object  motion  research  apply  directly  to  the  motion  of  fluids  and  solids. 

The  proposed  Image  Correlation  Velocimetry  (ICV)  method  that  will  be  de¬ 
veloped  in  the  present  discussion  has  roots  in  both  the  correlation  methods  used 
in  measuring  fluid  flow  and  the  deformation  of  solids,  outlined  in  Sec.  1.1,  and  the 
gradient  methods  used  in  measuring  optical  flows,  outlined  in  Sec.  1.2. 


1.1  Correlation  methods 

Several  techniques  for  determining  fluid  flow  velocities  from  particle  image  pairs 
{e.g.,  Willert  Sc  Gharib  1991)  employ  an  optimization  that  relies  on  some  form  of  a 
cross-correlation  function,  e.g.. 


max 

a 


j  Bo  (x)  £:,(«) 


(1) 


with 


^  =  x-l-a  , 


(2) 
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where  a  is  a  vector  parameter  to  be  determined  by  the  optimization  procedure  and 
A  is  tlie  correlation  region.  Tlie  distribution  of  the  image  irradiance,  E{x,t),  is 
known  at  times  to  and  ti,  i.e.. 


Eq{x)  =  E{x,tQ)  and  Ei(x)  =  E{xJi)  .  (3) 


The  average  velocity,  ,  within  the  correlation  region  is  then  approximated  by, 


u 


A 


a 

ti  —  to 


(4) 


The  drawback  of  this  method,  having  only  two  parameters  to  quantify  the  motion, 
is  that  it  cannot  resolve  displacements  properly  where  there  are  large  displacement 
gradients  within  the  correlation  region.  Anticipating  this  problem,  and  being  very 
interested  in  displacement  gradients,  rescarcliers  in  Solid  Mechanics,  apply  tech¬ 
niques  which  include  higher  order  deformations  of  the  displacement  field  within  a 
correlation  volume.  For  example,  McNeill  et  al.  (1987)  describe  a  method  whereby 
a  model  of  the  image  displacement  field  (mapping)  is  used  in  a  least  squares  opti¬ 
mization  procedure,  i.e., 


min  f  [  f;o{x)  -  £li(0  .  (5) 

a,Va 

The  affine  mapping, 

^  =  X -f  a -f  ( Va)  •  dx  ,  (6) 

is  used  as  an  example  of  such  a  function,  and  the  displacement  a  and  the  four 
components  of  Va  are  treated  as  parameters  to  be  determined  by  the  optimization 
procedure.  However,  any  physically  motivated  mapping  can  be  used  in  place  of 
Eq.6. 


In  both  these  methods,  the  image  data  are  integrated  over  a  region  and  require 
no  spatial  differentiation.  Since,  for  two-dimensional  images,  only  a  few  parameters 
are  extracted  from  the  optimization,  these  methods  are  relatively  immune  to  noise 
and  lend  themselves  to  fast  solutions. 


1.2  Gradient  methods 

A  method  for  determining  the  velocities  of  visual  features  in  an  image  was 
presented  by  Horn  &  Schunck  (1981).  This  visual  velocity  is  termed  “optical  flow” 
to  differentiate  it  from  the  velocities  of  (materiiU)  objects  in  an  image,  e.g.,  a  shadow 
moving  across  the  ground  has  a  perceived  velocity  that  is  markedly  different  from 
that  of  the  ground,  and  a  rotating  featureless  disk  will  have  no  visual  velocity  at 
all.  The  fundamental  eejuation  used  by  Horn  Schunch  to  determine  the  optical 
flow  was. 


^+u.VE  =  0, 


(7) 
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where  E  is  the  image  irradiance  (a  scalar)  and  u  is  the  optical  flow  velocity.  The 
diff’erential  terms,  dEfdt  and  VJ5,  can  be  estimated  from  the  image  data  and  the 
component  of  u  along  VE  \s  calculated  using  Eq.  7.  Methods  employing  equations 
of  this  type  are  called  gradient  schemes.  Note  that  no  velocity  can  be  calculated 
using  Eq.  7  if  there  are  no  features  or  gradients  in  the  image,  i.e.,  if  VE  =  0. 
In  addition,  because  Eq.  7  employs  only  the  component  of  u  along  VE.  velocity 
components  along  the  equi-scalar  contours  of  E  cannot  be  determined  using  Eq.  7 
alone.  This  limitation  was  designated  the  “aperture  problem”  by  Wallach  (1976). 
The  terminology  is  somewhat  misleading  and  is  only  used  here  in  reference  to  the 
convention.  More  appropriate  names  might  be  the  “characteristic  problem,”  or  the 
“direction  problem,”  because  the  problem  is  finding  the  velocity  along  the  charac¬ 
teristic  direction,  E  =  constant.  Gradient  schemes  also  have  the  problem  that  finite 
difference  approximations  of  the  spatial  and  tempered  derivatives  are  necessary.  A 
problem  with  such  approximations  for  the  derivatives  is  related  to  the  Nyquist  s£un- 
pling  criterion,  where  aliasing  in  the  image  data  can  effect  the  velocity  estimates. 
To  minimize  this  problem  in  taking  the  gradient,  the  motion  between  images  should 
be  less  than  half  the  smallest  local  spatial  scale,  A^;,  of  the  field,  i.e., 

|u|  ih-to)  ^  1 

a;  <  2 

(c/.  Eq.  7),  where  (tj  —  to)  is  the  time  between  images. 

The  uncertainty  of  the  so-called  aperture  problem  can  be  solved  in  some  cases 
by  applying  constraints  to  the  motion,  e.g.,  the  motion  is  of  a  rigid  body  (see 
Murray  &  Buxton  1990,  for  example),  or  a  limited  class  of  deformable  bodies  (see 
Terzopoulos  &  Metaxas  1991). 


Horn  &  Schunck  (1981)  applied  a  global  constraint  to  Eq.  7  (in  two  dimensions), 
by  solving  for  u(x,  t)  using  an  optimization,  t.e., 

f  i  ^  +  u  VE"!  -f-  )  d^x ,  (9) 

7a  \  J  J 

where  a  represents  the  constraint  cost  function  in  the  optimization  process,  cmd 
k  balances  the  relative  cost  of  a  and  Eq.  7.  In  particular,  Horn  &  Schunck  chose 
smoothness  £ts  a  constraint,  i.e.. 


The  idea  of  including  constraints  in  the  optimization  process  that  determines  the 
velocity  field,  over  an  area,  is  important  in  the  context  of  the  method  to  be  discussed 
below.  Note  that  the  constraint  in  Eq.  9  need  not  be  included  in  the  optimization 
integral.  Instead,  it  could  be  included  as  a  feature  of  the  optimization  technique. 
See  Murray  &  Buxton  (1990). 
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2.  Proposed  methodology 

A  succession  of  images  can  represent  anything  from  the  motion  of  cars  on  a 
highway  to  the  tremsport  of  a  dye  marker  in  water.  We  take  the  view  that  given 
successive  image  representations,  there  exists  a  transformation,  or  mapping,  of  the 
local  image  intensity  data  that  takes  one  image  to  the  next.  In  many  cases,  while 
the  equation  of  motion  of  the  imaged  field  may  be  known,  the  mapping  taking  one 
image  to  another  may  not  be.  Successive  images  combined  with  the  equations  of 
motion,  however,  often  allow  us  to  approximate  the  mapping. 

The  mapping  of  one  image  to  the  next  can  be  developed  by  considering  the 
Lagrangian  displacement  field  ^(x,t)  of  the  image  sequence.  Specifically,  if  x  is 
the  coordinate  of  a  point  on  an  image  at  some  initial  time  to,  ^(x, t)  represents 
the  coordinate  of  this  point  in  a  subsequent  image  recorded  at  a  later  time  t.  If  we 
imagine  the  image  sequence  as  the  result  of  a  continuous  recording  process,  we  can 
assign  a  Lagrangiiin  image-Row  velocity  field,  referred  to  as  “optical  flow”  in  the 
discussion  and  literature  cited  above,  t.e., 

u[^(x,t),t]  =  ,  (11) 

to  the  continuous  displacement  field  ^(x,t)  that  takes  an  initial  point  x  in  the 
image  recorded  at  time  to,  to  the  point  ^(x, t)  on  the  image  recorded  at  time  t. 
We  recognize  that,  for  the  case  where  the  images  represent  fluid  flow,  e.g.,  successive 
images  of  a  convected  scalar,  the  image-flow  velocity  field,  u(x,t),  may  be  quite 
different  from  the  fluid-Row  velocity  field  Uf(x,  t).  The  extent  to  which  the  former 
represents  a  good  approximation  for  the  latter  is  a  separate  issue  that  c^ln  only  be 
addressed  in  the  context  of  the  detmls  of  the  particular  imaging  process  and  the 
fluid-flow  field. 

In  the  proposed  implementation,  local  series  approximations  for  the  displace¬ 
ment  mapping  au’e  used  in  conjunction  with  an  integral  form  of  the  equations  of 
motion.  A  global  nonlinear  correlation  (optimization)  process  is  employed  to  esti¬ 
mate  the  image-flow  velocity,  vorticity,  deformation  rate,  etc.,  of  the  imaged  data 
field.  “Series,”  in  this  discussion,  will  denote  “Taylor  series.” 

In  the  context  of  fluid  mechanics  measurements,  we  will  focus  on  images  of 
continuous,  passive,  convected  scalars,  e.g.,  dye  markers,  carried  by  a  fluid.  As  will 
be  illustrated  using  a  pair  of  Voyager  2  images  of  Jupiter  (Sec.  5),  however,  any 
marker  in  the  flow  can  be  used. 

The  method  will  be  developed  for  three  dimensions  and  can  yield  three- 
dimensional  velocity  fields.  The  method  can  also  obtain  two-dimensional  velocity 
fields  from  images  of  two-dimension^ll  flows.  In  a  concession  to  the  limitations  of 
typical  data  acquisition  systems  today,  however,  the  method  will  be  applied  here 
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to  a  two-dimensional  flow  zind  also  to  two-dimensional  slices  of  three-dimensional 
flows.  We  note  that  in  some  cases,  two-dimensional  imaging  devices  can  be  used  to 
obtain  approximations  to  three-dimensional  image  data,  e.g.,  Dalim,  et  al.  1991.  A 
short  discussion  of  the  implications  of  correlating  two-dimensional  slices  of  three- 
dimensional  data  is  presented  in  Sec.  4. 


2.1  Fluid  displacement  and  equations  of  motion 

To  see  how  the  image-flow  velocity  field  can  be  calculated  from  three- 
dimensional  image  data  sets,  spaced  in  time,  first  consider  a  Lagrangicin  description 
of  a  flow  being  imaged.  Figure  1  illustrates  the  motion  of  fluid  particles  within  a 
volume,  V .  Fluid  elements  at  x ,  in  a  neighborhood  V ,  at  time  to ,  are  convected  to 
locations  ^(x,  t)  at  a  later  time  t.  The  displacement  field,  ^(x,  t),  can  be  thought 
of  as  a  transformation  of  the  field  x,  at  time  to,  to  ^{x, t).  Given  the  image 
displacement  field  ^(x,  t),  the  image-flow  velocity  field  is  then  given  by  Eq.  11. 


Fig.  1  Motion  of  a  fluid  volume. 


Using  this  Lagrangian  field,  ^(x,  t) ,  one  could,  in  principle,  integrate  the  equa¬ 
tion  of  motion  of  the  imaged  scalar,  t.e.. 


—  -I-  u  •  Vc  =  PV^c  , 
at 


(12) 


to  obtain 


rt\ 

ci[^(x,ti)] -co[^(x,to)] -P  /  V2c[^(x,t),t]  dr  =  0  , 

Jto 


(13) 
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where 

co(x)  =  c(x,<o)  and  ci(x)  =  c(x,ti)  (14) 

represent  the  c(x,t) -field  at  times  to  and  ti ,  respectively,  mid  V  is  the  appropriate 
diffusion  coefficient. 


In  the  first  two  examples,  the  motion  of  food  coloring  in  glycerine  (Sec.  3)  and 
dilute  fluorescein  in  water  (Sec.  4),  are  examined.  In  these  flows,  the  diffusion  of 
the  dye  mairkers,  in  the  time  interval  between  successive  images,  is  relatively  small 
and  may  be  neglected,  i.e.. 


P 


«:  1  , 


(15) 


where,  ti  —  to  =  0.1  s  is  the  time  between  images,  £  =  50  to  80 /zm  is  the  imaging 
resolution,  and  the  diffusion  coefficients  are  no  larger  than  V  =  10”®m^/s.  In 
addition,  we  note  that  the  Schmidt  number  is  large  in  both  flows,  t.e.. 


Sc  =  i^/V  >  10^  , 


(16) 


where  i/  is  the  kinematic  viscosity.  In  the  first  example  of  the  dye  marker  in  glyc¬ 
erine,  the  fluid  flow  is  two-dimensional,  as  is  the  image,  and  the  image-flow  velocity 
field,  u(x,  t),  is  a  good  representation  of  the  fluid-flow  velocity  field,  Uf(x,  t).  In 
the  second  example,  both  the  flow  velocity  and  the  imaged  scalar  field  are  three- 
dimensional,  while  the  image  is  two-dimensional.  As  we  will  discuss,  the  image-flow 
field  need  not  necessarily  represent  the  flow  velocity  field,  in  that  case.  In  the  third 
example,  the  motion  of  the  imaged  quantity  in  the  Jovian  atmosphere  (Sec.  5)  does 
not  follow  cmy  simple  equation  of  motion.  In  that  exaimple,  the  derived  image-flow 
velocity  field  can  be  expected  to  be  an  even  poorer  representation  of  the  fluid-flow 
velocity  field. 


In  cases  where  the  diffusion  of  the  imaged  scalar  cem  be  ignored,  Eq.  13  becomes 

ci[C(x,ti)]  -co[C(x,to)]  =  0  .  (17) 

Equation  17  represents  a  significant  simplification  over  Eq.  12,  its  differential  coun¬ 
terpart.  It  contains  no  spatial,  or  tempor£il,  derivatives  and  suffers  few  of  the 
drawbacks  associated  with  the  gradient  methods  discussed  eailier  (Sec.  1.2). 

Using  the  integral  equation  of  motion  (Eq.  17)  in  place  of  the  differential  equa¬ 
tion  of  motion  (Eq.  7)  in  the  optimization  (Eq.  9),  and  generailizing  the  optimization 
to  three  dimensions  then  yields  an  expression  for  determining  ^(x,  to)  and  ^(x,  ti), 
z.e.. 


mm 


d^X 


(18) 
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In  the  spirit  of  the  correlation  methods  discussed  in  Sec.  1.1,  where  the  type 
of  motion  within  the  correlation  volume  is  limited  to  translation  alone,  the  present 
method  restricts  ^(x,  to )  and  ^(x,  ti )  to  the  first  few  terms  of  a  series  approximation 
for  ^(x,  t) .  While  many  representations  of  the  displacement  field  could  be  employed, 
a  (Taylor)  series  representation  is  used  because  the  first  two  orders  in  the  series 
expansion  correspond  to  physical  fluid  mechanical  qumitities,  x.e.,  the  velocity  vector 
and  the  velocity  gradient  tensor.  More  importantly,  the  series  approximation  has 
the  additioned  benefit  of  enforcing  smoothness  in  the  displacement  and  displacement 
gradient  fields  within  a  correlation  volume. 


2.2  Displacement  field  and  kinematic  quantities 

In  the  case  of  fluid-flow  images,  the  quantity  ^(x,  <)  is  a  complicated  nonlinear 
function  of  the  imaging  process,  the  nonlinear  convection  dynamics,  and  x.  Local 
estimates  of  this  function  cein  be  made  by  Taylor  series,  expanding  ^(x,  t)  in  space, 
at  some  time  t,  in  an  image  correlation  volume,  V.  This  yields, 

^(x,0  =^(Xc;t)  +  (x-Xc)- V^(Xc;t)+  ^  [(x  -  Xc)  •  V]^  ^(Xc;  t)  +  . . .  (19) 

In  this  expression,  Xc  denotes  the  center  of  the  image  correlation  volume,  V,  at 
time  toi  and  V^(Xc;t)  denotes  the  gradient  of  ^(x,<)  with  respect  to  x,  evaluated 
at  Xc-  Figure  2  plots  the  number  of  parameters  used  in  the  optimization  process  as 
a  function  of  the  order  used  in  the  series  expansion,  for  two  and  three  dimensions. 
Figure  3  illustrates  the  eflFect  of  the  various  orders  of  the  expansion  on  a  two- 
dimensional  square  “volume.” 


Using  a  finite  difference  approximation  in  time  for  the  velocity,  Eq.  11, 2md  the 
series  representation,  Eq.  19,  evaluated  at  times  to  and  ti,  yields  an  estimate  for 
the  velocity  within  the  correlation  volume,  t.e., 


u[^(x;t),t]  ~ 


^(x;ti)  -^(x;  to) 
ti  —  to 


(20) 


where  to  <  t  <  t\.  Similarly,  taking  the  spatial  gr2idient  of  Eq.  20  yields  an  expres¬ 
sion  for  the  velocity  gradient  tensor  within  the  correlation  volume,  t.e.. 


Vu[^(x;t),t] 


[^(x;t),t] 


V$(x;ti)  -  VC(x; tp) 
t\  —  to 


(21) 


Vorticity,  divergence,  and  strain  rate  can  then  be  obtained  from  the  components  of 
the  estimated  velocity  gradient  tensor,  Eq.  21. 


Number  of  parameters 
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Order  of  series  expansion 


Fig.  2  This  plot  shows  the  rapid  increase  in  the  number  of  parameters  used  in  the 
optimization  procedure  with  increasing  order  of  the  scries  expansion.  “  *  ”  is 
for  a  2-D  expgmsion  and"©”  for  3-D. 

There  is  some  freedom  in  choosing  the  coordinate  transformation  at  the  initial 
time  to»  to)-  Dur  choice  is  to  have  the  coordinate  description  at  the  initial  time 
to  correspond  with  the  local  Eulerian  coordinates  at  that  time,  i.e., 

^{x,to)  =  x.  (22) 

In  terms  of  the  series  expansion,  Eq.  19,  this  means  that 


V$(Xc;to)  =  I 


(23) 


where  I  is  the  identity  tensor,  and  all  other  higher  order  derivative  terms  are  iden¬ 
tically  zero. 
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Translation  (no  deformation) 


Linear  deformation 


Quadratic  deformation 


Cubic  deformation 


Combined  translation  and  deformation 


Fig.  3  The  effect  if  translation  and  various  orders  of  deformation  on  a  two- 
dimensional  square  “volume.”  Sec  Sec.  2.2. 
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2.3  Seeking  a  global  solution 


A  solution  for  the  coefficients  of  the  series  expansion,  i.e.,  ^(Xc;t),  V^(Xc:t), 
etc.,  can  be  obtained  in  a  neighborhood  ai’ound  Xc  using  the  expansions  for 
^{x, to)  and  ^(x, t)  from  Sec. 2.2.  The  unknown  coefficients  of  the  series  ex¬ 
pansion,  ^(Xc;t),  V^(Xc;t),  etc.,  are  treated  as  parameters  in  an  optimiza¬ 
tion  process.  To  minimize  the  difference  between  two  data  sets  (images),  in 
a  least  squares  sense,  for  a  single  correlation  volume,  we  use  the  optimization, 
Eq.  18,  in  conjunction  with  the  series  approximations  developed  in  Sec. 2.2,  i.e., 

min  [  (  (  Cl  [£(xc;ti) -f  (x-Xc)  •  V^{xc;ti) -I- ...  ]  -  co[x]  I  -|- 

^(xo;ti).V^(Xc:«i) - tv  \  ^ 

-(-  ka"^  ^  d^x  .  (24) 

The  optimization  implied  in  Eq.  24,  combines  many  of  the  best  features  of  corre¬ 
lation  methods  and  gradient  methods,  while  eliminating  many  of  the  deficiencies. 
Specifically,  this  optimization  method  has  high  immunity  to  noise,  uses  the  equa¬ 
tions  of  motion,  can  incorporate  constraints,  requires  no  differentiation  to  calculate 
the  displacement  field,  and  can  capture  displacement  gradients  within  a  correlation 
volume. 


In  principle,  a  single  correlation  volume  covering  the  entire  image  and  a  series 
approximation  of  a  high  enough  order  can  be  used  to  capture  the  entire  image 
displacement  field.  In  practice,  however,  employing  a  series  approximation  beyond 
the  third  order  (cubic)  term  is  impractical  because  of  the  rapid  increase  in  the 
number  of  parameters  in  the  optimization  process  with  increasing  order  (see  Fig.  2), 
and  the  associated  increase  in  the  computationeil  time  and  complexity.  In  the 
present  calculations,  when  the  quadratic  term  is  not  sufficient  to  capture  the  image 
deformation  over  the  entire  flow  field  using  a  single  volume,  as  is  usually  the  case, 
several  series  expansions  residing  in  smaller,  adjacent,  correlation  volumes  are  used 
in  place  of  the  single  large  volume. 


To  construct  a  global  optimization  using  a  number  of  local  series  expansions,  we 
require  that  neighboring  correlation  volumes,  with  independent  series  expansions, 
must  yield  consistent  results.  In  the  present  method,  we  use  the  expansion  for  the 
displacement  field  about  one  correlation  volume  to  estimate  those  of  its  neighbors. 
The  displacement  field  of  these  neighbors  is  also  estimated  in  terms  of  their  own  se¬ 
ries  expansions.  The  root-mean-square  difference  between  displacements  estimated 
by  neighboring  correlation  volumes  is  applied  as  a  constraint  cost  function.  Since 
it  is  necessau'y  to  refer  to  a  number  of  scries  expansions,  it  is  useful  to  define  the 
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/  4,iX42  ; 


^,(*42 ;  ‘) 


krO 


Fig.  4  The  fcl  -  A;8  solid  circles  denote  the  points  in  V,  used  by  the  constraint  cost 
function  <7^ .  The  empty  circles  denote  their  counterparts  estimated  by  the 
neighbors  Vji  -  Vj4 . 

taylor  series  for  ^(x,  t)  in  a  neighborhood  centered  about  Xc, ,  as 

^,(x;  t)  =  C(Xc. ;  t )  +  (x  -  Xc. )  •  V^(xc. ;  0  +  ^  [(x  -  Xc. )  •  V]^  ^(Xc. ;  t )  + . . .  .  (25) 

When  series  expansions  about  multiple  points  are  then  employed,  the  minimization, 
Eq.  24,  is  modified,  i.e., 

min  f  (  {ci[^,{x;ti)J  -  co[x]  j  +  kaf  )  d^x  ,  Vi  .  (26) 

This  minimization  is  performed  within  all  the  Vi,  simultaneously,  and  the  square 
of  the  constraint  cost  function, 

|^.(xfc;0 -C>{xfc;0  r  (27) 
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is  applied  to  provide  global  continuity  of  the  solution,  denotes  the  series  expul¬ 
sion  about  the  “j”  neighbors  of  V^,  and  X/t  denotes  the  “A:”  points  of  comparison 
between  the  solutions  and  See  Fig.  4.  In  the  present  method,  eight  points 
about  each  correlation  volume  axe  used  for  comparison,  three  with  each  of  four 
neighbors.  These  eight  points  are  sufficient  to  define  the  series  expansion  coeffi¬ 
cients  of  the  correlation  volume,  up  to  quadratic  order. 

The  present  implementation  of  the  method  can  solve  Eq.  2G  for  two-dimensional 
flow  up  to  the  cubic  term  in  the  local  series  expansions,  but  the  series  is  usually 
truncated  at  quadratic  order.  The  optimization  of  Eq.  24  is  accomplished  using  a 
multidimensional  minimization  process,  with  image  data  between  pixels  estimated 
using  bilinear  interpolation.  See  for  example  Press,  et  al.  (1988). 


2.4  Minimization  parameters  in  two  dimensions 


Typical  CCD  imaging  technologies  today  are  limited  to  two-dimensioncil  (spa¬ 
tial)  data.  This  is  not  a  problem  if  the  flow  being  imaged  is  also  two-dimensional. 
This  section  describes  how  the  method  is  applied  in  two  dimensions.  First,  we  de¬ 
velop  the  terms  of  the  series  expansion,  Eq.  19,  for  two-dimensional  flow.  With  the 
two-dimensional  vector 

=  x-Xc  ,  (28) 

as  the  position,  x ,  relative  to  the  center  of  the  correlation  volume,  Xc ,  the  terms 
of  the  series  expansion  at  a  time  ti  appear  as  a  constant  term. 


:itl  1  —  , 

L«oJ 


where  the  are  the  vector  coordinates  of  the  center  of  the  correlation  volume  at 
the  time  ti ,  t.e., 

«i  =  ^i(Xc;ti)  ,  (30) 

a  linear  term, 

(x-Xc)  •  VC(Xc;ti)  =  ,  (31) 

Q2.1  0:2,2]  [02 

where  the  Oij  represent  the  first  order  deformations  and  rotations  of  the  image 
field  within  the  correlation  volume,  t.e.. 


O-i  j  — 


_  c)^,(Xc;ti) 


a  quadratic  term. 


2!  [0:2,11  0:2,12 


0:1,22 

02,22 
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where 


a  cubic  term. 


.  =  L  (  ^  ^ 

^i,jk  —  21  yj  ^  j.  _  2  J  dxjdxk 


(34) 


-[(x-Xc)  •  vf  $(xc;ti)  = 


with 


«1,111  ai.222 

«2.111  ^2,112  02.122  02.222 


r  1 

J2<52 

6xSl 

L  51  J 

(35) 


.  =]l(  3 

3!  \i+j  +  k-3)  dxj 


dxi^dxi 


(36) 


and  so  on  for  higher  order  terms.  The  otijk  ^^nd  ocijki  are,  respectively,  related 
to  the  second  and  third  derivatives  of  the  displacement  field  within  the  correlation 
volume,  i.e.,  by  Eqs.  34  and  36. 

The  velocity  and  velocity  gradient  (Eqs.  20  and  21),  can  also  be  written  in 
terms  of  the  parameters  of  Eqs.  29  and  31  and  the  series  expansions  at  times  to 
and  ti ,  i.e.. 


1 

Ql  -  Xc 

LvJ  ti  —  to 

.  ^2  -  2/c  . 

where 


and 


Xc  = 


Xc 

.Vc] 


Vu(Xc,<)  = 


du/dx 

du/dy 

1 

Ol.l  -  1 

Ol,2 

dv/dx 

dv/dy 

ti  —  to 

02,1 

02,2  —  1  . 

(37) 


(38) 


(39) 


Alternatively,  the  velocity  gradient  can  be  written  in  terms  of  the  in-plane  vorticity 
and  rate-of-strain  tensor,  i.e., 


Vu(Xc,0  = 
where  ujz  is  the  vorticity,  i.e.. 


0  — 

a;,/2  0 


SlI  Sxy 

Sxy  Syy  J 


(40) 


u;,  = 


02,1  ~  Ol.2 

t\  —  to 


(41) 


2md  Sij ,  Syy,  and  Sxy,  are  the  components  of  the  rate-of-strain  tensor,  i.e.. 


Ql.l  -  1 
tl  —  to 


(42) 
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«yy  = 


Q2.2  —  1 

tl  —  tQ 


and 


Sxy  -  2 


1  0:1,2  +  02,1 


i\  —  to 


(43) 


(44) 


An  interesting  quantity  to  consider  is  the  second  invariant  of  the  rate-of-strain 
tensor  (see  Cantwell  1992,  for  example),  i.e.. 


2  +  2  s^y  +  s 


yy 


(45) 


3.  Couette  flow  between  concentric  cylinders 

An  apparatus  to  generate  a  Couette  flow  between  concentric  cylinders  was 
fabricated  for  the  purpose  of  testing  the  method.  The  cylinders  were  made  from 
248  mm  lengths  of  stock  Plexiglas  tubing.  The  inner  and  outer  radii  of  the  annular 
region  between  the  cylinders  were  nominally  25.2  mm  and  40.9  mm.  The  cylinders 
were  from  stock  Plexiglas  tubing,  so  the  uncertainties  in  the  radii  were  ±lmm.  The 
outer  cylinder  was  rotated  with  a  rotation  rate  of  1.1  rad/s,  with  the  inner  cylinder 
stationary.  In  this  example,  employing  a  dye  marker  in  glycerine,  the  fluid  flow  is 
nominally  two-dimensional  and  the  marker  follows  the  flow.  Hence,  the  image-flow 
velocity  field,  u,  can  be  accepted  as  a  good  representation  of  the  fluid-flow  velocity 
field,  Uf. 

Images  were  recorded  using  a  Texas  Instruments  Multicam  MC-1134GN  Multi- 
Mode  B/W  Camera.  The  data  were  stored  digitally  using  an  in-house  multiple  frame 
grabber  (12-bit  A/D),  designed  by  Dan  Lang  and  Paul  Dimotakis  of  GALCIT,  set  to 
record  up  to  28  of  the  1134  x  468  pixel  gray  level  images  from  the  camera,  spziced 
by  100  msec  (sidjustable  between  33  and  267  msec).  Because  the  horizontcd  and 
vertical  spacing  of  the  pixels  were  not  equal  on  this  CCD,  grid  spacings  and  image 
correlation  volumes  with  a  ratio  of  1:1.74  (vertical:horizont2il  pixels)  were  used  to 
yield  a  uniform  sp^lcing  of  the  data  in  the  read  image  plane.  Flow  visualization 
WEIS  performed  by  randomly  distributing  red  food  coloring  (dye)  on  the  surfaice  of 
the  fluid.  To  provide  bsicklighting  for  the  dye  mEirker,  the  fluid  beneath  the  surfsice 
contEiined  a  trsuislucent  white  suspension  of  3  /xm  aluminum  oxide  ( AI2O3  )  ptu-ticles 
in  glycerine.  When  illuminated  from  the  side,  this  provided  nesirly  uniform  white 
bsw^klighting  for  the  dye  being  imaged  on  the  surface.  Because  of  the  depth  of  field 
of  the  imaging  and  the  high  density  and  uniform  distribution  of  the  Eiluminum  oxide, 
scattering  from  individual  pEirticles  was  not  detectable  in  the  video  images. 


Fig.  5  Initial  |)laccnicnt  of  scries  expansion  neighborhoods.  Each  square  denotes  a 
control  volume.  The  small  circle  at  the  center  of  each  control  volume  denotes 
the  center,  or  control  point,  of  a  scries  expansion. 


In  the  present  investigations,  only  the  outer  cylinder  was  rotated,  hence  the 
velocity  field  can  be  written  as  (c.g.,  Schlichting  1979), 


Ur{rJ)  _  ,  uoir,9)  rfr^-r^lr 

— -  =  0  and  — - =  — - - j-  , 


(46) 


where  r ,  6,  Ur.  and  uo  are  the  radial  and  angular  positions  and  velocities,  respec¬ 
tively,  is  the  rotation  rate  of  the  outer  cylinder,  and  r,  and  Tq  are  the  inner 
and  outer  radii  of  the  cylinders.  In  this  flow,  the  divergence  is  zero,  i.c.. 


V  u  =  0  , 


(47) 


and  the  vorticity  is  uniform,  i.c.. 


2110 

1  -  (rifr„f 


uJz(r.O)  =  ej  Vxu(r.ff)  = 


(48) 


I  ’K (I  I  )is|)liii  «‘iiiriil  t»(  serifs  iifi'^liIxMluMMls.  lOOins  l:ili*r.  alhtwiii*',  only  I  raiislat  ion. 

I'ur  (liis  siini>le  lost  (  asc.  nilif  rorn'la! ion  volnincs  spacfii  hy  If)  jiixols  vntirally 
anti  TS.i}  jiixfls  lioii/onl ally,  weif  used  To  ea]»tnre  the  How.  d'lie  lofatioiis  nf  the 
iiiia;’,'-  fttiTflal  ion  vtthinies  at  the  initial  time  are  dei)ieted  in  I'i^.  T).  d'he  results 
of  the  t  (uri’lat ion  jiroecss.  allowiiif^  <»iily  displacfincnt  of  tin*  forrolation  volume 
(/,ert)th  order  sei  ies  ex|»ansion)  is  shown  in  Fi^- h.  Fij;ure  7  tlemonst  rates  liow  the 
How  is  ea|)ture<l  wImmi  the  correlation  proee.ss  is  extcmletl  to  inelmh*  hijHiei'  oriler 
terms  in  the  series  expansittn.  d'o  tjuantifv  this  improvement,  the  value  of  the 
mini  mi/at  ion  funet  ional.  not  including'  the  eoiit  rihut  ion  of  t  he  const  raint  s.  is  plot  t  ei| 
in  I'i^;.  ;ts  a  function  of  the  order  of  the  series  expansion.  As  can  he  .seen,  the 
f;realesl  improvement  is  reali/.ed  with  the  introduction  of  linear  deformations  in  the 
etirrelal  ion  process. 

The  results  iisiiif’,  the  lorrelation  method  are  compared  with  the  t  heuret  ical 
( analy  t  i<  ill )  two-dimensional  values  for  the  vortieity  ami  diverj^enee  in  I'ahle  1. 
Sim  e  holh  the  vortieity  and  the  diveif^t’iiee  are  uniform  in  the  analyti<  al  solution, 
only  a  sin;’.le  value  is  presented.  'The  uncertainty  in  the  theoretical  vortieity  is  a 
result  of  the  eeeent  l  icit  v  of  t  ln'  cylinders  used  in  the  experiment .  'I'lie  iincei  taint  ic> 
in  the  experimental  values  are  one  standard  deviation.  Hecaiise  this  ijow  h.u-  a 
nearly  lineai  velm  ity  jirolile.  we  saw  only  small  i  hani'.es  in  these  estimale-.  hevonil 
tlie  linear  order.  Some  of  the  “uni'ert  aiiity '  in  tin'  experimental  result.-  for  the 


Fig.  7  Displacement  of  series  neighborhoods  estimated  using  higher  order  terms  in 
the  correlation  process. 

vorticity  reflect  the  expected  variations  of  the  vorticity  within  the  flow. 


Theoretical: 


Experimental: 


a;,  [sec  *] 


V  •  u[sec 


-3.6  ±0.1 


0.0 


-3.6  ±0.3 


-0.1  ±0.1 


Table  1.  Comparison  of  theoretical  and  experimental  vorticity  and  divergence. 
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Order  of  series  expansion 


Fig.  8  Difference  between  Couette  flow  images  under  the  mapping,  quantified  by 
the  value  of  the  optimization  functional  (arbitrary  units),  as  a  function  of 
the  number  of  terms  in  the  two-dimensional  series  expansion. 
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4.  Cylinder  wake  flow.  Two-dimensional  slices  of  three-dimensional  data 

Ill  this  section,  we  present  the  results  of  applying  the  two-dimensional  cor¬ 
relations  to  two-dimensional  slices  of  a  three-dimensional  flow  in  the  wake  of  an 
impulsively  started  circular  cylinder.  Here,  the  cylinder  is  1.75  cm  in  diameter  and 
45.5  cm  long.  It  is  drawn  through  a  distribution  of  fluorescein  dye  in  water  at  a 
speed  of  1.27  cm/s.  The  Reynolds  number  in  this  case  is, 

Ud 

Re= - %  220  ,  (49) 

1/ 

where  U  is  the  cylinder  speed,  d  is  the  cylinder  diameter,  and  v  is  the  kinematic 
viscosity. 


Fig.  9  Initial  placement  of  series  expansion  neighborhoods.  Each  square  denotes  a 
control  volume.  The  small  circle  at  the  center  of  each  control  volume  denotes 
the  center,  or  control  point,  of  a  series  expansion. 


The  CCD  camera  and  data  acquisition  system  are  the  same  ais  for  the  Cou- 
ette  flow  test  case  (Sec.  3).  Laser  sheet  illumination  is  provided  by  a  Continuum 
model  YG661-10  frequency-doubled  YAG  laser.  The  laser  was  operated  at  532  nm, 
300  mJ,  5  ns  pulses  width,  at  a  rate  of  10  Hz.  The  flow  here  is  three-dimensional  in 
both  the  velocity  field  zmd  scalar  distribution. 


Fig.  10  Displacement  of  grid  after  100  ms,  estimated  using  the  nonlinear  correlation 
process. 


Figures  9-14  demonstrate  the  method  on  images  of  a  vortical  structure  forming 
in  the  wake  of  the  cylinder.  These  images  were  tadcen  after  the  cylinder  had  traveled 
about  8  diameters.  The  image  at  the  initial  time  is  shown  in  Fig.  9,  and  100  ms  later, 
in  Fig.  10.  In  this  case,  the  series  approximation  used  in  the  correlation  process  Wcis 
expanded  to  quadratic  order.  Figure  11  shows  the  displacement  of  the  centers  of 
the  correlation  volumes. 


The  two-dimensional  vorticity  is  displayed  in  Fig.  12.  A  large  vortical  region 
can  be  seen  in  the  wake  of  the  cylinder.  The  two-dimensional  divergence,  presented 
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Fig.  11  Displacement  of  centers  of  grid  over  100ms. 


in  Fig.  13,  exposes  the  three-dimensionality  in  the  flow.  Figure  14  plots  the  second 
invariant  of  the  rate-of-strain  tensor.  Note  the  region  of  strain  (rate)  that  seems 
to  follow  the  periphery  of  the  large  vortical  structure.  This  could  be  a  region  of 
vorticity,  from  the  previously  shed  vortical  structure,  that  is  being  strained  around 
the  current  one. 


As  a  general  observation,  an  important  issue  arises  when  imaging  a  two- 
dimensional  (planar  image)  slice  of  a  three-dimensional  field  of  a  continuous  scalar, 
c(x,t),  as  in  the  previous  example.  An  out-of-plane  component  of  the  fluid-flow  ve¬ 
locity,  Uf ,  coupled  with  an  out-of-plane  component  of  the  scalar  gradient,  Vc(x,  t) , 
will  contribute  to  the  in-plane  image-flow  velocity  u.  In  this  case,  the  equation  for 
the  in-plane  image  flow  can  still  be  written  as. 


dc  dc  dc  dc  dc 


=  0 


(50) 


where  we  have  assumed  that  the  image  irradimice  E(x,t)  is  proportional  to  the  two- 
dimensional  slice  of  the  scalar  concentration.  c(x,  t),  and  where  the  in-plane  image- 
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Fig.  12  Contours  of  constant  plane-normal  vorticity,  ^  x  u  =  Solid 

contours  denote  zero-vorticity.  Long  dashes  denote  positive  values  and 
short  dashes  negative.  Contours  spaced  by  0.5 s"^ . 

flow  velocity,  u  =  {u,v)  —  d^/dt  (Eq.  11),  is  the  one  derived  from  the  minimization 
function,  as  described  above. 


Considering  the  transport  of  the  three-dimensional  iso-scalar  surfaces  (see 
Fig.  15),  we  find  tliat  the  two-dimensional  u  =  (u,v)  in-plane  image-flow  velocity 
components  are  related  to  the  three-dimensional  Uf  =  {u{,Vf,W{)  fluid-flow  velocity 
and  the  three-dimensional  scalar  gradient  components.  In  particular,  we  have, 

_  dc  dc/dx 

^  ^  ^  ^  dz  {dc/dx)'^  +  {dc/dy)^ 

dc  dcjdy 

Uf  -f  W(  [dcldx)"^  d-  (defOy)'^ 


As  can  be  seen,  by  substituting  Eq.  51  in  Eq.  50.  these  relations  recover  the  three- 
dimensional  transport  eciuation  for  a  conserved  scalar  field  c(x.t),  in  the  ciisc  of 


I*'lG.  Contours  of  constant  divergence,  V  -  u  =  §7  + 
<*ontours  ar<^  spacnul  by  0.5 s~* , 


As  in  Fig.  12.  tlje 


negligibh^  diffusion,  i.c.. 


Or. 

Ot 


+  Ilf  • 


Oc 

Ox 


Oc  Oc  Oc  Oc 


=  0  . 


These  re.sults  provide  us  with  tlu;  criteria  for  when  t!u'  in-plane  iniage-llow 
v<'locity  (  an  be  r(?gard(Ml  as  a  good  approximation  to  the  in-plaiu'  fluid-flow  vt>locit y. 
In  particular,  we  will  have. 


and 


r  Or/O, 

lllf  ' 

[  iOc/Ox)-  +  («>c/(b/)-  , 

-»l 

Hlf  ' 

<  Oc/Oz  > 

\  ihi 

^  uh-joj-r-  +  / 

I  On  ■ 

(52a) 
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Fig.  14  Contours  of  constant  second  invariant  of  the  rate-of-strain  tensor  (Eq.  45). 
=  [(§7)^  +  -  I  ]  •  Contours  spaced  by  0.5 s'^. 


We  can  see  that  if  W{  is  small,  or  if  dcjdz  is  small,  or  both,  by  the  measure  in  Eq.  52, 
then  the  in-plane  image-flow  velocity  field  can  be  accepted  as  a  good  representation 
of  the  in-plane  fluid-flow  velocity  field. 

Finally,  since  this  method  estimates  the  in-plane  image-flow  velocity  u,  and 
not  the  fluid-flow  velocity  Uf,  its  application  to  two-dimensional  image  slices  of 
three-dimensional  scalar  field  data  is  identical  to  its  application  to  two-dimensional 
scalar  data. 
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Fig.  15  Two-dimensional  slices  of  a  three-dimensional  scalar  field  c{x,t).  Uf  in¬ 
dicates  the  3-D  fluid-flow  velocity  and  u  is  the  resulting  two-dimensioned 
in-plane  image-flow  velocity. 

5.  Voyager  2  images  of  Jupiter 

The  method  is  also  illustrated  on  a  pair  of  the  images  of  the  atmospheric  dy¬ 
namics  of  Jupiter  teiken  by  Voyager  2.  These  images  were  taken  from  the  “Voyager 
Time-Lapse,  Cylindrical- Projection  Jupiter  Mosziics,”  by  Avis  &  Collins  (1983). 
640  X  350  pixel  subimages  of  rotations  349  and  350  were  used  in  the  correlation 
process.  The  subimage  spans  168°  to  97°  longitude  and  0°  to  —46°  latitude  (the 
equator  is  at  the  top  of  the  image).  The  subimage  from  rotation  349  is  shown  in 
Fig.  16,  with  an  overlay  of  the  initial  placement  of  the  correlation  volume  neighbor¬ 
hoods.  The  vertical  line  on  the  left  is  a  reference  line  which  is  to  be  deformed  using 
the  mean  zonal  velocities  of  Jupiter  from  Limaye  (1985).  Figure  17  shows  the  same 
region,  one  rotation  later,  with  the  associated  grid  deformed  by  the  nonlinear  cor¬ 
relation  method.  On  the  left  is  the  reference  line  from  Fig.  16,  carried  by  the  mean 
zonal  flow.  The  displacement  of  the  centers  of  the  correlation  volumes  is  shown  in 
Fig.  18. 
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Fig.  16  Initial  placonient  of  correlation  volnines  ovcrlayecl  on  a  snb-iinage  of  ro¬ 
tation  349  from  the  "Voyager  Time-Lapse.  Cylindrical-Projection  Jupiter 
Mosaics.”  Each  scpiare  of  the  grid  denotes  a  correlation  volume.  The  ver¬ 
tical  line  on  the  left  is  a  reference  line  to  be  carried  with  the  mean  zonal 
velocity  of  Jupiter  (see  Fig.  17). 


Fig.  17  Deformation  of  the  correlation  volumes  (see  Fig.  16).  after  one  rotation. 
The  line  on  the  left  wiis  deformed  from  the  vertical  line  in  Fig.  16  using  the 
mean  zonal  velocity  of  Jupiter. 
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Fig.  18  Displacement  of  grid  control  points  (centers  of  correlation  volumes)  after 
one  rotation.  The  lines  on  the  left  denote  the  displacement  via  the  mean 
zonal  flow  of  Jupiter. 


6.  Conclusions 

Series  expansions  of  image  displacement,  in  conjunction  with  a  global  nonlinear 
correlation  method,  can  be  used  to  measure  fluid  velocities,  and  velocity  gradients, 
from  pairs  of  continuous,  convected,  scalar  images.  It  is  shown  that  increasing 
the  order  of  the  expansion  can  improve  the  accuracy  of  the  results.  The  proposed 
method  does  not  require  discrete  particles  and  may  also  be  used  in  situations  where 
there  is  a  natural  marker  alrecwJy  in  the  flow,  e.g,,  species  concentration  can  be 
used  to  measure  velocities  in  compressible  flows.  The  method  is  developed  for 
three-dimensional  data  sets  and  demonstrated  on  two- aimensional  images  of  fluid 
flow. 
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